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Abstract: Using combined strong coupling and hopping parameter expansions, we derive 
an effective three-dimensional theory from thermal lattice QCD with heavy Wilson quarks. 
The theory depends on traced Polyakov loops only and correctly reflects the centre sym- 
metry of the pure gauge sector as well as its breaking by finite mass quarks. It is valid 
up to certain orders in the lattice gauge coupling and hopping parameter, which can be 
systematically improved. To its current order it is controlled for lattices up to Nr ~ 6 
at finite temperature. For nonzero quark chemical potentials, the effective theory has a 
fermionic sign problem which is mild enough to carry out simulations up to large chemical 
potentials. Moreover, by going to a flux representation of the partition function, the sign 
problem can be solved. As an application, we determine the deconfinement transition and 
its critical end point as a function of quark mass and all chemical potentials. 
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1 Introduction 

The determination of the QCD phase diagram is a fully non-perturbative problem, because 
QCD is strongly coupled on the scales relevant to heavy ion collisions and astrophysics, 
i.e. for temperatures T < 400 MeV and baryon chemical potentials fiB ~ — 3 GeV. On the 
other hand, a direct first principles approach by Monte Carlo simulations of lattice QCD 
is ruled out by the sign problem, with a complex fermion determinant for quark chemical 
potential = /xb/3 7^ prohibiting importance sampling. Existing workarounds based 
on reweighting, Taylor expansions in fi/T or simulations at imaginary chemical potential 
followed by analytic continuation all introduce additional systematic errors and require 
fi/T<l in order to be valid. For an elementary introduction, see [1]. As a consequence, 
the QCD phase diagram remains largely unknown. 

This situation warrants additional investigations by effective theory methods. Popular 
approaches are based on models which share the chiral and/or the Z(3)-symmetry with 
QCD, such as PNJL models, Polyakov loop + quark meson models, sigma models etc. For 
recent discussions and references, see [2-4]. Other approaches start from QCD directly 
and employ Dyson-Schwinger [5] or functional renormalisation group methods [6], using 
particular truncations. A general difficulty with effective theories is to assess the associated 
systematic errors. 

Recently, a systematic derivation of a 3d centre-symmetric effective theory for finite 
temperature SU{Nc) Yang-Mills by means of a strong coupling expansion has been pre- 
sented, followed by numerical simulations [7]. The effective theory depends on Polyakov 
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Figure 1. Integration of thie cube consisting of 6 plaquettes. In the middle, four spatial link 
variables have been integrated over and we end up with two doubly occupied plaquettes. Integrating 
two of the remaining Ui yields the structure on the right which gives a factor [Tr(J7'''i7)]^ = iV^. 

loop variables only, and can be improved order by order. Its couplings Aj = Xi^ji^N^) 
are calculable functions of the lattice gauge coupling /3 and the temporal lattice extent 
Nt of the original finite temperature lattice theory. The influence of the various couplings 
can be checked. Remarkably, the effective theory with only one coupling reproduces the 
correct order of the deconfinement transition for S'C/(2), ^[/(S), and moreover permits a 
quantitative estimate of the critical couplings f3c{Nr) for sufficiently fine lattices, such that 
a continuum extrapolation of Tc is feasible. 

In this work we extend the effective theory to include heavy but dynamical fermions 
of mass M by means of a hopping parameter expansion. This theory permits to explore 
the phase diagram of QCD with heavy quarks in the (M, T, //) parameter space. A similar 
approach to the fermionic sector was taken in [8], which however left the gauge sector in 
the original 4d form. Here we extend the effective fermionic action to order in the 
hopping expansion. As we shall see, for /x = we once again obtain good agreement with 
full 4d simulations where such results exist. However, the 3d effective theory allows for a 
solution to the sign problem and to explore the full range of quark chemical potentials with 
numerical ease, making contact with the region of asymptotically large chemical potentials. 
As a first application of the theory, we map out the entire deconfinement critical surface 
delimiting the region of first-order deconfinement transitions in the (Mu^di Mg, jj) parameter 
space. 

In Section 2 we derive the effective theory and discuss its numerical evaluation, compar- 
ing a direct simulation with Polyakov loop degrees of freedom with that of a flux represen- 
tation free of the sign problem. Section 3 is devoted to the study of the QCD deconfinement 
transition, followed by our conclusions and a discussion of further prospects in Section 4. 

2 The effective action 

2.1 Finite temperature SU{Nc) Yang-Mills theory 

For the paper to be self-contained and to fix the notation, we briefly summarise the main 
formulae for the SU{Nc) pure gauge case [7]. Starting point is the partition function for 
the 4d Euclidean Yang-Mills theory with Wilson gauge action. 



Finite temperature and the bosonic degrees of freedom imply the use of periodic boundary 
conditions in the compactified time direction with A'^ slices, setting the temperature scale 
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by the lattice spacing a, T = l/(aA^T-). An effective three-dimensional theory emerges by 
integrating out the spatial degrees of freedom to get schematically 



[dUo]exp[-Scs] , 



-5, 



eff 



In J [dUi] exp ^ Yl Up + Tr f/j) 



XiSi + A25'2 + 



(2.21 



where the effective couplings A, are functions of the original parameters, Aj = Xi{/3,Nr). 
The introduction of the logarithm is convenient in order to employ the graphical cluster 
expansion, described e.g. in [9], which features only connected diagrams. 

An important observation is that only graphs winding around the time direction are 
needed for the effective action. To see this, consider for example the integration of a cube 
shown in Fig. 1, which is a valid graph with non- vanishing contribution. Nevertheless, since 
it does not wind around the temporal lattice extent, the spatial link integrations suffice 
to remove all dependence on link variables. As a result, the cube contributes only as a 
function to the effective partition function and hence cancels in expectation values or 
renormalised quantities such as the physical free energy density. This is true for all graphs 
which do not wind around the lattice. 

After spatial integration, the interaction terms Si in Eq. (2.2) then depend on the link 
variables only via Polyakov loops 



L(x) = Tr W{x) = Tr Yl Uo{x, r) . 



(2.3) 



T = l 



Thus we transform the remaining path integration measure from temporal link variables to 
traced Polyakov loops, which introduces a reduced Haar measure denoted by . We now 
consider SU (3) and parametrise the measure by two angles of the diagonalised Polyakov 
loop, providing another factor e^, so that at every spatial lattice site x we have 



He, 



dW-- 



y = I In ( 27 - 18lLp + 8Re(L3) - 



T = l 



dLe^ 



de 



d(j) e 



2V 



(2.4) 



In Eq. (2.2), we arrange the interaction terms in ascending order of their leading power 
of /3 in the strong coupling expansion of the corresponding effective coupling. The first such 
interaction term is between nearest neighbours {ij) in the fundamental representation and 
has the form 



(2.5) 



where Lj = L{xi). This leading order contribution comes from N-^ temporal plaquettes 
building a chain around the lattice, followed by spatial link integration. Knowledge of 
the relations Xi{f3,Nr) allows to convert the critical couplings Aj^c of the three-dimensional 
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Figure 2. Critical temperature for the pure gauge deconfinenient transition, extracted from tlie 
3d effective theory, Eqs. (2.7, 2.6). 



theory to those of the original theory. In particular, for the nearest neighbour coupling 
Xi{u, N-r) we find 



Ai(u,2) 



exp 



4 

u exp 



Xi{u,Nr > 5) 



u exp 
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(2.6) 



where the character expansion coefficient u = u{(3) = /3/18 + . . . of the fundamental 
representation is used instead of /3 due to better apparent convergence. 

Altogether Eq. (2.2) has infinitely many couplings with loops in all irreducible repre- 
sentations, to all powers and at all distances. In the strong coupling region, the Si are 
the more suppressed the higher the index i. In [7] we saw indeed that the influence of the 
next-to-nearest neighbour coupling as well as the adjoint coupling are negligible within the 
current level of accuracy when investigating the phase transition. We thus neglect these 
and higher-order correction terms in the following. Summing up powers of the nearest- 
neighbour interaction term [7], the effective theory for thermal Yang- Mills reads 



^eff 



\ i / <ij> 



(2.7) 
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Figure 3. Different terms occurring in the hopping expansion. Left: Plaquette. Middle: An 
example of a 6-link graph. Right: Generalised Polyakov loop, i.e. the loop winds around the 
temporal direction n times before the link variables are traced. 

The deconfinement transition of Yang-Mills theory is reflected in an order-disorder transi- 
tion of the effective theory. We determine the critical coupling to be Ai^c = 0.18805(2) = Aq, 
cf. Sec. 2. 5^, and, by inverting Eq. (2.6), extract the fSd^r) which agree with full 4d Monte 
Carlo results within better than 10% accuracy up to Nr = 16. Using the non-perturbative 
beta-function for a(/3) provided in [10], these can then be converted to deconfinement tem- 
peratures, as shown in Fig. 2. Note that all points stem from a single determination of X\ c 
in the effective theory. A continuum extrapolation in a^, i.e. is feasible and predicts 

a deconfinement transition of Tc = 250(14) (50) MeV. The second, systematic error is taken 
as the difference between the 0{u^^) and 0{u^) in Eq. (2.6). Encouraged by this result, 
we now proceed to extend the effective theory to include heavy Wilson fermions. 

2.2 Heavy fermions: LO hopping parameter expansion 

Heavy fermions are conveniently introduced using the hopping parameter expansion, which 
is described in [9] at zero temperature and discussed in [11, 12] for finite temperature. The 
quark part of the action for Nf mass-degenerate flavours with masses Mf = M can then 
be written as 

-S, = -Nfy'^TrH[U]\ ^1 = ^7 , (2.8) 

1=1 

with the hopping matrix 

H[U]y,^ = Y,^y,x+i>{i + lu)Uy{x) , = . (2.9) 

±u 

Thus each hop to a neighbouring lattice site gives a power of the hopping parameter k. 
The quark chemical potential n is introduced as usual by a factor e"'^ (e~"'^) multiplying 
link variables in positive (negative) time direction [13]. The Kronecker delta in Eq. (2.9) 
requires that the graphs in the hopping expansion be closed. In Fig. 3 we show several 
graphs appearing in the hopping expansion. As an example and in order to establish a 
physical mass as reference scale, we calculate the pion mass to leading orders, 

aM^ = -21n(2K) - 6k^ - 54^"^ - 2Ak^ + 0(k^ , k^u^) . (2.10) 

1 — u 

At finite temperature there are also graphs with a nontrivial winding number, such as 
generalised Polyakov loops 

TrH^"(f) = Tr rQc/o(x,r) , 1 < n < oo , (2.11) 

^The slight discrepancy between Ao and the Ai,c of [7] is due to finite-size effects, as in the latter 
determination system sizes only up to A'^s = 14 were employed. 
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winding around the lattice n times before being traced, cf. Fig. 3 (right). 

We obtain the effective action from the full partition function in the same way as in 
pure gauge theory, 

Z = J [dUo][dU^] exp [Sg -Sg] = J [dUo] exp [Ses] , 

-S,s = lnj [dU^] exp [Sg - S^] . (2.12) 

We are now faced with a double series expansion in u{f3) and k, i.e. the effective couplings 
depend on both parameters and A'^^. Furthermore, quarks of finite mass lead to terms in 
the action which explicitly break the Z{3) symmetry present in the pure gauge case. The 
effective action may then be written as 

oo oo 

-S,s = Y, Ai(^x, K, Nr)St - 2Nf k, /x, N^)S^ + h,{u, k, ft, N^)Sf'^] . (2.13) 



2 = 1 



i=l 



The Aj are defined as the effective couplings of the Z(3)-symmetric terms Sf, whereas the 
hi multiply the asymmetric terms Sf. Consequently, only the latter are /^-dependent and 
we recover pure gauge theory for k ^ 0, as in the full theory. We have not included the 
factor 2Nf in the definition of the hi, since there are Nf mass-degenerate quarks with 2 
spin degrees of freedom, all giving the same contribution. The hi and hi are related via 



hi{u, K, fl, Nr) = hi{u,K,-fl,Nr) . 



(2.14) 



We shall now derive combined strong coupling and hopping parameter expansions of 
these effective couplings, again by employing the graphical cluster expansion. Similar to 
the case of pure gauge theory, graphs contributing to the cluster expansion have to wind 
around the lattice in the compact r-direction. Hence, the leading order contributions are 
Polyakov loops, and we can read off the leading-order couplings hi and hi: 

hiSl + hiSl'^ = -Yl [{'iKe''^)^^L{x) + (2Ke-"^)^-L*(x)] ^ hi = (2Ke'^'^)^^(2.15) 

Note the minus sign which is due to the anti-periodic boundary conditions for fermions. It 
is possible to sum up the contributions of all generalised Polyakov loops oriented in positive 
time direction, 



expi-2iV^^^ 



X n=l 



(-ir 



(2Ke"'^)"^-Tr {W 



n 



1 + (2Ke' 



aii\NT 



w 



2Nf 



,(2.16) 



using exp Tr In A = detA, and similarly for the conjugate loop. The three-dimensional 
effective action to leading order in the hopping expansion then corresponds to the static 
approximation and reads 



Zeff = J [dUo] n + 2AiReL*L,] ^|^JJ det 



l + hiWs)il + hiW, 



2Nf 



.(2.17) 
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Figure 4. Left: Graphs reducing to Polyakov loops after spatial link integration resulting in 
0{k^v}) corrections, with 1 < / < Nr — 1. Right: Graphs of C'(k^^^+^), leading to interactions 
between Polyakov loops at distance a after spatial link integration. 



2.3 Heavy fermions: beyond leading order 

Corrections to the static approximation come from graphs that also contain spatial link 
variables. To order 0{k'^) there is the plaquette term, see Fig. 3 (left). Its contribution 
can be absorbed in the gauge term leading to a shift in j3 and hence to a K-dependence of 



u{l3) n(/3, k) . 



(2.18) 



Higher-order graphs containing 6 and more links [9] , for which Fig. 3 (middle) is an example, 
may be neglected to the orders to which we have calculated our effective couplings. 

Next, we consider leading-order corrections to the winding graphs in the full (3+1)- 
dimensional action and observe their effect on the 3d effective theory, cf. Fig. 4. Graphs 
on the left lead to 0{k'^) corrections to the couplings hi, hi after spatial link integration. 
We have calculated the first correction to be 



Nr-l 



{2Ke'"'f-6NrK^ ^' = (2Ke"^)^^6A^,K' 



,u — u 



1-u 



(2.19) 



by summing over 1 < Z < N-r — 1 possible attached plaquettes. As a result we get the 
higher order version of the effective coupling, 



hi = {2Ke 



aii\NT 



1 + GNrK- 



,u — u 



Nr 



+ ... 



l-u 

Including all corrections up to 0{u^''k"^), with n + m = 7, we obtain 

l-^z^^-i 



hi{u,K,Nr > 3) = {2Ke^^)^- 



exp 



6NrK^U 



1 - u 



Ak 



(2.20) 



(2.21) 



where we used the fact that terms ~ or higher can be resummed by writing the 
correction as an exponential. The graph in the lower right of Fig. 4 contributes a k- 
dependent term to the already included nearest-neighbour interaction of Eq. (2.5). This 
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exemplifies how Xi{u, Nr) becomes \i{u, k, A'^^). Note that these corrections do not change 
the form of the effective action Eq. (2.17). 

By contrast, a graph that gives rise to a new term in the effective action is shown 
in the upper right corner of Fig. 4. Together with its oppositely oriented counterpart it 
implies interactions between Polyakov loops of the same orientation, 

h2S^ + hS^'^ = (h2LiLj + h2L*L*') , (2.22) 



<^3> 

h2 = (2^e»^)2^^ ^ , h2 = (2/^e-'^^)2^^ ^^"^^ 



However, the terms arising from both graphs on the right of Fig. 4 are parametrically of 
high order O (k^^'^"''^) and will thus be neglected along with even higher-order corrections. 

For the remainder of this work we thus approximate thermal lattice QCD with heavy 
quarks by the effective theory Eq. (2.17) with the couplings hi of Eq. (2.21) and Ai of 
Eq. (2.6) with u{/3,k) as in Eq. (2.18). Since we retain only one coupling of each sort, let 
us drop the index 1 and use the index i only to label lattice sites to lighten the notation. 
Finally, we rewrite Z^s in a form more convenient for its numerical evaluation. A useful 
identity for the determinants is 



det(l + hW) = 1 + hTrW + h^TvW^ + = 1 + hL + h^L* + h 
det(l + hW^) = 1 + hTiW^ + h'^TvW + h^ = l + hL* + J?L + h 



The heavy quark contribution per flavour and site xi then becomes 

Qi{h,h) = [(1 + hLi + h^L* + h^){l + hL* + h^U + h^)]^ , (2.23) 
and the 3d effective theory for thermal QCD with heavy quarks reads simply 



eflf 



/"( J]dL, e^>Qf^ (/i,/i) J J] (l + 2AReL,L*) . (2.24) 

\ i / <ij> 



In our numerical investigations we consider the partition function Eq. (2.24) for Nj = 1. 
Since the deconfinement transition at high temperature happens at small /i, we can recover 
an arbitrary number of flavours by using the approximation 

det(l + hW)^i ^ exp (NfhL) , h{Nf) « —h{Nf = 1) . (2.25) 

Nf 

2.4 Observables 

In our numerical simulations we are interested in the phase structure. With a standard 
action linear in its couplings, one would typically construct observables using energy- and 
magnetisation-type fields defined as 

Eiin = ^ 2ReLjL* , Qnn = -^| (2.26) 

* <ij> i 
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on finite volumes V = N^. In our case the non-linear action implied by Eq. (2.24) suggests 
to use a different definition, closer to what actually drives the dynamics of the system, 

* <ij> ^ i 

These are proportional to the previous definitions in the limit of small couplings {X,h,h). 
For the finite size scaling analyses of the phase structure we used the susceptibility and the 
Binder cumulant constructed from the observables O € {E,Q, Enn^Qim}, 

Xo = NU{0')-{Of), B,,o = llo-{0)Py ■ ^^-^^^ 
2.5 Flux representation and numerical simulation 

The system described by the partition function Eqs. (2.17, 2.24) has a complex action 
and thus suffers from a sign problem for h ^ h {/j, > 0). In the simplified case where 
the local variables Li take only values in the center, Li E Z(3), i.e. the Potts model in a 
magnetic field, the sign problem can be solved by using a cluster algorithm [14] or a change 
of variables to obtain a flux representation [14, 15]. The latter approach can be generalised 
to the case of SU (3)-valued Polyakov loops as was done in [16] for a related model. The 
flux representation of the partition function reads 

ZesiX,h,h)= ^ Jl Wb{nb,mb)'[[Wi{pi,qi,ni,mi) . (2.29) 

As the 5C/(3)-valued variables have been integrated out, the degrees of freedom are now 
local incoming and outgoing currents, nj^^,mj^^ = 0,1, respectively, located on the links 
connecting the site i and i + fi, fi = ±1, . . . , id as well as charges (or monomers) ni,mi = 
0, . . . , 4 located on the sites of our cubic lattice. Wb{nb, rub) is given by 



Wb{nb,inb) 



1, if (nb,mb) = (0,0) 

A, if (nb,mb) = (0,1) or (1,0) (2.30) 
otherwise . 



*m 



If we rewrite the factor Q in Eq. (2.23) as a power series in L,L*, Q = Yin ■m^n,mL"'L 
then the site weight Wi reads 

W,{pi,qi, n,,m,) = ^n.,m, J dU L'?^+-^L*f»+"^^ > , (2.31) 

where the last expression contains a traced SU{3) link integral given in closed form in [17, 
18]. Wi is positive for all // > and the model no longer has a minus sign problem. Locally, 
this integration creates a Potts constraint Qi + rii = pi + nii mod 3, where 

qi= ^ rii^^ , Pi= ^ mi,fi ■ (2.32) 

fi=±l,...,±d /i=±l,...,±d 

Thus, local currents are conserved modulo 3. 
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Figure 5. Z^g/Zes as a function of A for the pure gauge case, crossing at Aq « 0.18805. 




The partition function Eq. (2.29) is well-suited for an application of the worm algo- 
rithm [19] and its variants. To enable sampling of configurations in the presence of an 
external magnetic field {h, h 0) we implemented a variant of the algorithm presented 
in [20]. Having changed the degrees of freedom, the observables Eqs. (2.27) have to be 
re-expressed. The Polyakov loop and its complex conjugate are now given by 

= ^|^ln^- W/^), :^^E^l)^ W/^), (2.33) 

where (n) and (m) denote the average number of monomers of either type. The relations 
become exact only in the limit h,h ^ 0, due to the non-exponential form of the quark part 
of expression Eq. (2.17). The quark density is given by 

The worm algorithm relies on the sampling of the 2-pt function G{i,j) = {LiL*) rather 
than the partition function Eq. (2.29). G{i,j) can be estimated during a worm update [19] 
as G{i,j) = h{i,j)/Z, where h{i,j) corresponds to the histogram of configurations with 
the worm head (say L*) at site j and its tail (L) at site i. 

An observable more suitable for the flux representation is the free energy of an inter- 
face enforced in the broken phase by twisted boundary conditions in the e.g. z-direction, 
Li+Nss^ = e*'^Lj with (p = ^q,q = 0, 1,2. The partition function with twisted boundary 
conditions is 

\nx,mxt 

where jz = Ylx&P iP'zix) — mz{x)) is the flux through the plane P = {x \ z = Ns — 1}. If 
we consider the case h = h = 0, i.e. the model representing pure gauge theory, then the 
spontaneous breaking of centre symmetry is signalled by the ratio 

Z,s I 0, for T > Te . ^ ^ 
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Figure 6. Left: Average sign, Eq. (2.37), of the fermion determinant in the effective theory in the 
vicinity of the critical point for various chemical potentials on a 2A^ lattice. Right: Quark density 
calculated with Zcs from Eq. (2.24) (Metropohs) or Eq. (2.29) (worm) on 24^ lattices. 



In Fig. 5 we show Z^/Zeg as a function of A for (/) = 2tt/3. Since our model has a weak 
first order transition several volumes cross at the transition point Aq = 0.18805(2). 

On the other hand, compared to the full theory the sign problem in the standard 
representation of the effective theory, Eq. (2.24), is mild, similar to the case of the Potts 
model in an external field [21]. Using the modulus of the determinant in a Metropolis 
algorithm while reweighting in its phase, system sizes up to Ng = 24 for values of the 
chemical potential up to fi/T ~ 3 can be reached, with a fully controlled average sign. To 
demonstrate this, we monitor the ratio of the full and phase quenched partition functions, 

(sign)|| = = e-^'^-^^''") , Zeffi I : phase quenched. (2.37) 

^effll 

The corresponding difference in free energy density is a volume-independent measure for 
the strength of the sign problem and is shown in Fig. 6 (left). Even for the largest system 
sizes to be used in our scaling analyses, the average sign remains significant and fully 
controlled up to large chemical potentials /u/T ~ 3. This is corroborated by comparing a 
physical observable such as the quark density between the worm algorithm (without sign 
problem) and standard Metropolis algorithm with Polyakov loops and reweighting in the 
phase of the determinants. Fig. 6 (right) shows that no difference is discernible between 
the two ways of evaluating the observable. 

We thus conclude that the sign problem can be fully controlled and even solved for 
our effective theory. Since the observables of interest are more readily accessible in the 
original degrees of freedom, we have mainly used the Metropolis algorithm for the following 
numerical simulations. 



3 The deconfinement transition for heavy quarks 
3.1 Zero baryon density 

As a first application of the effective theory we investigate the deconfinement transition of 
QCD with heavy quarks as a function of quark mass and chemical potential. We begin by 
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Figure 7. Schematic phase diagram for ^ = 0. Left: Nature of the QCD thermal transition for 
different quark masses. Right: Expected phase diagram for the effective theory corresponding to 
11 = 0. 

considering the case of zero baryon density, shown schematically in Fig. 7 (left). In the 
pure gauge limit, the deconfinement transition is of first order. Dynamical quarks at any 
fixed A'^^ break the global Z(3) symmetry of the QCD action explicitly. As a consequence, 
the phase transition weakens with decreasing quark masses until it vanishes at a critical 
point. For still lighter quarks the deconfinement transition is an analytic crossover. 

This behaviour is inherited by the effective theory. For a given Nf and fi = 0, we have 
h = h and the effective theory has two couplings, (A,/i). The first-order phase transition 
of the one-coupling theory extends to a first-order line with a weakening transition as h 
increases. Eventually the transition should vanish at a critical point, as sketched in Fig. 7 
(right). These expectations are based on the known results of the 3d 3-state Potts model 
in an external field [14, 21], which shows the same symmetry breaking pattern. While the 
behaviour of the system in the vicinity of the critical point is dictated by the universality 
class, the location of the transition in parameter space, in particular the critical parameters 
where it changes its nature, are not. Hence, our investigation will give valuable additional 
information on QCD. Previous investigations to locate the critical heavy lattice quark mass 
have been made on coarse A',- = 4 lattices for Nf = 1 [22], and in [23] for several flavours. 

In order to determine the phase diagram Fig. 7 (right), we follow a two-step procedure. 
First we determine the phase boundary, i.e. the pseudo-critical line Xpc{h) in the two- 
coupling space of the effective theory. In a second step, using dedicated finite size scaling 
analyses, we determine the order of the transition along that line, and in particular the 
location {Xc,hc) of the critical point. In order to accomplish the first task we fix the 
external field variable to the values h = 0.0002 — 0.0012 on lattice sizes A^^ = 10 — 24, and 
then scan for the corresponding pseudo-critical coupling Ape- As indicators for the phase 
boundary we use maxima of susceptibilities and minima of Binder cumulants constructed 
from the observables given in Eqs. (2.26, 2.27). We extrapolate these to infinite volume 
using 

Xp,{h, Ns) = Xpcih) + ci{h)/N^ . (3.1) 
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Figure 8. Pseudo-critical couplings ioi h = h = 0.0006. Left: B^ q^.^ on various volumes, the 
minima approach Xpc{h) in the thermodynamic limit. Right: Extrapolation of Xpc{h, N^) from 
several observables to infinite volume. 
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Figure 9. The pseudo-critical line or phase boundary in the thermodynamic limit. Also the critical 
point, Eq. (3.5), is shown. 



For each coupling and system size we generated at least 10^ configurations. Fig. 8 shows 
the minima of the Binder cumulant B^^q^.^ for h = 0.0006 and the extrapolation of these 
values along with those of other observables. This results in the pseudo-critical line shown 
in Fig. 9. It is well described by a hnear fit, due to the small magnitude of h and the 
argument given in [14]: along the line of first order transitions the free energy densities 
of the disordered (confined) phase and ordered (deconfined) phase are equal, fc{\,h) = 
fd{X,h). Expanding both sides about the pure gauge transition, {\o,h = 0), and noting 
that dhfc = {L) = in zero external field, we obtain 

dhfd 



^pc{h) = Ao - aih, ai 



dxifc - fd) 



(Ao.O) 



(3.2) 



A fit with = 0-26 yields ai = 1.797(18) and Aq = 0.18805(1) in very good agreement 
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xq 


Xq 
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0.000 My(lj 


0.000 1 L[L ) 


0.000 iZ[Z) 


V 


0.63(1) 


0.630 (fixed) 


0.64(2) 


0.630 (fixed) 




2.00(1) 


1.998(1) 






k 


16.88(4) 


16.87(2) 


1.58(1) 


1.58(2) 


fl 


-162(12) 


-157(3) 


7.6(8) 


7.0(2) 


f2 


460(60) 


(fixed) 


1(1) 


(fixed) 



Table 1. Critical end point he from fits to the scaling forms Eqs. (3.3) using data from Ng — 
20,22,24 lattices (all with acceptable x^/dof < 1.5). The known 3d Ising critical exponents are 
= 1.962(3) and v = 0.6302(1) [24]. 

with the estimate from the flux representation, see Fig. 5 and Sec. 2.5. An alternative 
approach to pin down the pseudo-critical line Xpe{h) would be to determine the principal 
axes £, A4 of the joint probability distribution P{E, Q) of our variables. The critical line 
is then defined by a certain symmetry condition on P{A4), demanding the vanishing of the 
third moment {A4^) = 0. We have checked explicitly around the critical point that the two 
approaches give consistent results. 

In order to locate the critical end point of the first-order line and to establish its 
universality class, we study finite size scaling of the data taken along the pseudo-critical 
line. Close to criticality our observables should scale according to 

XQ = N]/''f^^{x), B,,Q = fB,Jx), x^{h-he)Nl/'' (3.3) 

with dimensionless scaling functions /o, provided we move along the tangent hpc{X). In 
the vicinity of the critical point the scaling functions may thus be expanded, 

fo{x) = fo + fioo + f2X^ + ---, (3.4) 

which is the form to which we fit our data. We simulated lattice sizes Ng = 20, 22, 24 with 
statistics of ~ 7 • 10^ configurations per parameter set and used binning analyses to control 
autocorrelations. For the susceptibility XQi find 7/1/ consistent with the expected 3d 
Ising value, see Table 1 and Fig. 10. The fit was repeated fixing f2 = 0, u = 0.630 (3d 
Ising), and varying the fit range |x| < 0.01,0.02,0.03, with an overall stable outcome. For 
the Binder cumulant -64,(5, /o should approach a value characterising its universality class, 
(/o = 1.604 for 3d Ising [24]). The same fitting procedure as above was then applied to the 
-84,(3 data for Ns > 20, with compatible values for he and /o, i', cf. Table 1. The collapse of 
the data onto a universal curve under the appropriate rescaling is also shown in Figs. 10, 
11. 

We are then ready to identify the critical point in Fig. 9 (right) 

(^Xe = 0.18672(7), he = 0.000731(40)) . (3.5) 

The values of Xc,he can be converted into those of the couplings (3c, using Eqs. (2.6, 
2.21). In order to compare with previous work, we approximate Me/T with the relation. 
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Figure 10. Left: XQ as a function of hpc{X) for several volumes. Right: Same, now with rescaled 
axes, Eq. (3.3), to produce a collapsing curve described by the universal scaling function (see Table 
1). The vertical line marks the critical h, the horizontal line is /q. 
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Figure 11. Left: B^ q as a function of hpc{X) for several volumes. The dashed horizontal line 
denotes the 3d Ising value 1.604 which is the approximate crossing value of the largest volumes 
Ns = 20, 22, 24, the vertical line marks the critical point. Right: Same, now with rescaled horizontal 
axis. The horizontal line marks B^^q^O) as found numerically. Also shown is the universal scaling 
function obtained by the fit listed in Table 1, last column. 



Nf 


MJT 


^c{Nr = 4) 


Kc(4), Ref. [23] 


Ate (4), Ref. [22] 


1 


7.22(5) 


0.0822(11) 


0.0783(4) 


~ 0.08 


2 


7.91(5) 


0.0691( 9) 


0.0658(3) 




3 


8.32(5) 


0.0625( 9) 


0.0595(3) 





Table 2. Location of the critical point for /i = and Nr — 4. The first two columns report our 
results (we used for consistency the leading-order relation Eq. 2.15), the last two compare with 
existing literature. 
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valid for heavy fermions to leading order in the hopping expansion [25] 



/ M^ 
exp ( - -) 



h 



(3.6) 




The results are collected in Table 2 and are in reasonable agreement with the corresponding 
ones from simulations of 4d QCD with Wilson fermions [22, 23] at A^,- = 4. 

As in the case of pure gauge theory, our mapping of the critical effective couplings to 
those of QCD can in principle be done for any A'^,-, thus providing predictions for larger A',- 
which have not yet been simulated in 4d. However, before doing so we need to check how 
far we can trust our hopping parameter expansion. Fig. 12 (left) shows the predictions 
of the effective theory for Kc{Nr) to the orders k^, k^, resummed and unresummed. Also 
shown is the chiral critical hopping parameter, defined by the vanishing of the pion mass 
Eq. (2.10), and evaluated for the critical gauge couplings, Kch {uc{Nt-)) ■ Since we are 
expanding around infinite quark masses, self-consistency requires Kc ^ i^ch- Whereas the 
leading order soon crosses the Kch hue, the corrections are significantly below the leading 
order and the exponentiated versions further improve on this. Furthermore, literature tells 
us that Kch{P = 0) = 0.25 [26] and K,ch{P — ^ oo) = 0.125 [9], when all orders are taken into 
account. Hence, we conclude that Nr = 6 is the finest thermal lattice for which our Kc is 
still significantly smaller than Kch evaluated at the same gauge coupling, and has not yet 
crossed the continuum-extrapolated Kch- This is corroborated by the pion mass Mt^{uc, Hc) 
evaluated at the critical point, shown in Fig. 12 (right). We observe that beyond Nr ^ 6 
the differences between the non-trivial orders and grow larger, indicating that we 
leave the convergence region. 

These findings are to be contrasted with /3c(A^t) of the pure gauge effective theory, 
which are within a 10% range from the known 4d results up to A',- = 16. However, this 
is quite natural, since we have only three non-trivial orders in the hopping expansion, 
which are additionally truncated at a low order in n, compared to five orders in the strong 
coupling expansion. While we hope to extend our results to larger A'^,- by going to higher 
orders in k, at the moment we cannot take a continuum limit in the fermionic sector but 
consider our results valid up to A',- ~ 6. 

Within this range of validity, we may now discuss the sensitivity of the deconfinement 
critical line in Fig. 7 (left) on the cut-off. The critical pion mass in units of temperature 
marking the boundary of the first order deconfinement region shrinks slightly from A',- = 4 
to A't- = 6. This effect is even smaller in absolute units (Tc also decreases with increasing 
A',-), in contrast to the critical pion mass evaluated on the chiral critical line, which shrinks 
by almost a factor of two [27]. Higher orders in the hopping expansion are needed for a 
definite statement in our case. 

3.2 Finite baryon density 

We now study the deconfinement transition at finite baryon density. For fi ^ 0, we have 
h ^ h and need to consider the full parameter space of the effective theory, {\,h,h). 
The diagram in Fig. 7 (right) turns three-dimensional, with a surface of first order phase 
transitions terminating in a critical line. Since we are interested in the change of the critical 
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Figure 12. Critical hopping parameter Kc{Nr) (left) and critical pion mass MTr^c{Kc,u)/T (right) 
given by Eq. (2.10) for the end point of the Nf = 1 effective theory to different orders of the hopping 
expansion. 



quark mass with chemical potential, we prefer to map out the critical line by fixing different 
chemical potentials and then scan for the critical k. It is thus convenient to introduce the 
parameter 

h = /le"^/^ (= (2k)^" to leading order in k) (3.7) 

and to present our data in the parameter space (A, h, jJ-fT). 

For our simulations, we used values of fi/T = 0.1, . . . ,3.0 on three lattice sizes Ng = 
16, 20, 24. Data were produced at a given set (A, h, n/T) close to the critical point and later 
reweighted to near-by values of the couplings. Over 700 k configurations were produced 
for each parameter set and lattice size. For each chemical potential, the pseudo-critical 
line \pc{h) was identified as the curve of local minima in -64,^, Fig. 13, showing again 
linear behaviour as already observed for = 0. Indeed we can repeat the steps leading 
to Eq. (3.2), this time setting fc{X,h, n/T) = fd{X,h, n/T) along the first order line. For 
small fields h we expand about the Yang-Mills limit {\,h, fi/T) = (Ao,0,0) and obtain 



Xpc(h,fi/T) = Xo + ai{fj,/T)h, ai 



d\{fc - fd) 



cosh(^/r). (3.8) 

(Ao,0,0) 



Fitting to this form (x^/dof ~ 1.5) gives an intercept Aq = 0.18802(2) consistent with Aq 
in Eq. (3.2) and a slope 

ai{fi/T) = Ccosh(^/T) (3.9) 

with C = —1.814(3), which is compatible with Eq. (3.2). 

As for /X = 0, the critical points XdhcifJ-fT), fi/T) can be found again by evaluating 
B4^Q along each //-line on different volumes. To avoid doing the entire finite size analysis 
for all parameter sets we take the critical point as the crossing of the Ng = 24 data with 
the theoretical value of 1.604. The difference between this procedure and the crossing of 
individual volumes serves as an estimate for finite size effects. On the resulting critical 




3 - 3.5 4 0.0002 0.0004 0.0006 0.0008 0.001 

h X10-' h 



Figure 13. Left: Binder cumulant B4^Q{h,X) for chemical potential n/T = 1.4 {Ng = 24). The 
thick straight lines mark the locus of the minima at Ns = 16, 24, the surrounding contour lines 
correspond to the values 1.58,1.64,1.74. Right: Transition lines \pc{h, n/T) for several values of 
n/T between (top) and 3.0 (bottom). 



line, \c{hc{^^/T), ^/T) shows only weak dependence on /i/T within our statistical accuracy 
and varies around Ac ~ 0.18670(5). Note that this remains true even for large /i/T, 
as we shall see in Eq. (3.12).^ We exploit this behaviour to find a simple parametric 
description of the critical line in terms of the parameters of the original QCD action. 
Setting AA = Ac(/i/r) — Aq ~ const, we may rewrite Eq. (3.8) at the critical point for fixed 
/x/T as 

Un/T) = = • (3.10) 

Gcosh(/i/r) cosh(;u/r) 

A fit of all /z > data to Eq. (3.10) with only D as parameter performs indeed very 
well, yielding D = 0.00075(1) with x^/dof = 0.6. The corresponding hdfJ- = 0) = D 
is also compatible with our findings at /i = 0. The resulting curve is shown in Fig. 14 
(right). While the data seem to be reasonably well described by the Ansatz Eq. (3.10), 
a systematic underestimation at small fi/T hints towards a more complex law. Indeed 
our curve is the result of a first order expansion in {\,h,h), Eq. (3.8), followed by the 
approximation AA ~ const and a fit over the whole //-range. 

On the other hand, asymptotically large chemical potentials in the original lattice QCD 
are described by the limit 

K ^ 0,/i — > oo with Ke^^'^ = const . (3.11) 

The corresponding effective theory has two parameters, (A, h,h = 0). The critical point in 
this case is easily found by the same techniques, 

(Ac, /ic)|/i=o = (0.18668(2), 0.00142(2)) . (3.12) 



^Interestingly, the same observation is made in the 3d Potts model, where the spin coupling as a function 
of the external fields is nearly constant along the critical line, e.g. [21]. 
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Figure 14. Left: B4 Q for different system sizes along the pseudo-critical line for n/T = 1.4. The 
horizontal lines mark the thermodynamic limit value of B4 Q at criticality and the actual value of 
the crossing. Right: Critical line hc{fJ-/T); the curves are the fit to the cosh~^(/i/T) behaviour, 
Eq. (3.10), and the large-/z asymptotic limit, Eq. (3.11). The separate /i = determination is 
compatible with these findings. 



Using the leading order expression for h, Eq. (2.15), this gives the h = 0-critical curve 
hc{fJ-/T) which is also plotted in Fig. 14. Already for fi/T>1.5 the data are accurately 
described by the asymptotic density limit. Taking this result together with the curve 
Eq. (3.10) we thus have obtained a description of the Nj = 1 deconfinement critical line 
for all real chemical potentials! 

3.3 Imaginary chemical potential 

The QCD phase transitions and its limiting critical surfaces possess an analytic continu- 
ation to negative or imaginary chemical potentials, /j, = ifii. This has been exploited 
previously [28, 29], since in this case the fermion determinant is real positive and prop- 
erties of the critical surfaces can be calculated without sign problem. In particular, the 
deconfinement critical surface is expected to terminate in a tricritical line at /ij/T = 7r/3 
[30]. This value of imaginary chemical potential marks the boundary to an adjacent Z(3) 
centre-sector of the partition function [31]. Our effective theory correctly reflects the centre 
symmetry and its breaking in QCD, and hence the related phase structure. In this section 
we explicitly compute the continuation of the critical quark masses, i.e. the deconfinement 
critical surface, from /i = to n/T = in/S. 

Now also the fermionic part of our effective theory, Qx, is explicitly real. Numerically 
we follow the same approach as for real choosing values of Hi/T = 0.1 — ^, followed 
by determinations of the pseudo-critical and critical couplings. We observe increasing 
numerical difficulties as /ij approaches the boundary to the next centre sector. Moving 
along the critical line towards the Roberge- Weiss tricritical point, a crossover between 3d 
Ising and tricritical scaling sets in, thus obscuring the finite size analysis and demanding 
ever larger volumes. Controlled errors were obtained up to fii/T <0.8. 
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Figure 15. Left: Critical coupling he for imaginary chemical potential ^i, also shown is the analytic 
continuation of the real-/Lt fit to Eq. (3.10). Right: M^/T for Nf = 1 at both imaginary and real 
chemical potential. The curves represent Eq. (3.10) (and its analytical continuation), the large-/^ 
asymptote from Eq. (3.11), and the tricritical scaling, Eq. (3.13). 



With increasing ^Uj, the endpoint of the corresponding first order line is shifted towards 
higher hdni/T). The resulting critical line is shown in Fig. 15 (left). The pseudo-critical 
lines \pc{h, fii/T) develop a curvature for increasing /ij and thus invalidate a first order 
expansion of the free energy as done in Eq. (3.8). As a consequence, an analytic continu- 
ation of the real-/x fit Eq. (3.10) to imaginary chemical potential leads to a less satisfying 
description of the data. However, we may put real and imaginary /i data together and plot 
Mc/T as in Fig. 15 (right). (To convert h M, we use Eq. (3.6), with hce~^^'^ instead of 
simply he). It was demonstrated in [30] for the Potts model and strong coupling QCD that 
the critical quark mass at imaginary chemical potential is governed by tricritical scaling, 
with a scaling region extending all the way to real jJL. We thus attempt the corresponding 
two-parameter fit to tricritical scaling. 



T \T^) T 



3 J +It 



(3.13) 



which is shown in Fig. 15 (right). Fitting solely the region fi^ < yields K = 1.55(3), and 
Mtric/T = 5.56(3) with Xrcd ~ 0-13. Different numbers of flavours can be re-introduced as 
for = 0, obtaining 

= {5.56(3), 6.25(3), 6.66(3)} , for Nf = {1, 2, 3} . (3.14) 

Remarkably, the scaling function correctly describes the data up to large real chemical 
potentials, in fact well into the region described by asymptotically large /x. 

3.4 The deconfinement transition for all parameter values 

Collecting our results from the previous sections, we can describe the entire deconfinement 
critical surface from imaginary chemical potentials all the way to the large- /i limit. More- 
over, we have a simple and accurate parametrisation of the surface by stitching together 
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Figure 16. The deconfinement critical surface of QCD with heavy quarks. 

tricritical scaling, the cosh~^(/x/T)-behaviour for moderate /u and the curve describing the 
asymptotic Umit. Converting to a (2 + l)-flavour setting with different masses M^^d / Mg, 
we may plot the critical surface for the upper-right corner of the Columbia plot (Fig. 7, 
left). The final result is shown in Fig. 16. To convert to quark mass M/T we used Eq. (3.6), 
valid to leading order in the hopping expansion and for small N^- < 6. 

By fixing the quark content, e.g. Nf = 2, we can similarly draw the full deconfinement 
transition in (^, ^)-space, with Tq the pure gauge transition temperature. Since our 
quarks are very heavy, they give negligible contribution to the beta-function and we use 
again the Sommer parametrisation from pure gauge theory [10] to set the scale. As for 
11 = 0, our results are only controlled up to Nj- ~ 6. This yields the phase diagrams in 
Fig. 17. The deconfinement temperature for quarks of large but finite mass is almost fi- 
independent and smoothly approaches a constant value as expected for the quenched limit 
of infinitely heavy quarks. 

4 Conclusions 

We have applied strong coupling and hopping parameter expansions to lattice QCD with 
Wilson fermions at finite temperature and quark chemical potential. The resulting three- 
dimensional effective theory depends solely on traced Polyakov loops, i.e. complex scalars, 
with its couplings given analytically in terms of the original parameters of the theory, 
PjKjNt- This leads to a considerable simplification of numerical simulations. In the pure 
gauge limit we know five non-trivial orders in the strong coupling expansion, and the 
resulting critical parameters are valid up to Nj- ~ 16, enabling a continuum extrapolation 
of Tc- However, the hopping expansion has only been performed up to order k^, such that 
our theory with quarks is only robust to Nj- ~ 6 so far. Nevertheless, this corresponds 
to finer lattices than have been directly simulated in the heavy quark regime at the time 
of writing, and moreover offers intriguing features absent in the full 4d theory. The sign 
problem is mild enough to directly simulate the model, and can be solved completely 
within a flux representation of the partition function. Hence a numerical study for all 
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Figure 17. Left: Phase structure of the theory in ^)-space, for Nj = 2 and Nj- = 6. 

Above the surface the theory is deconfined. The critical hne (dashed) separates the cross-over 
(hght) and the first-order surface (dark). Right: Phase diagram of the theory for AI^/{2T) ~ 8.9, 
Nf = 2 and Nr = 6. 

values of the chemical potential is feasible. We have demonstrated this by computing 
the entire deconfinement critical surface. While this region of the QCD parameter space 
is far from the physical parameter values, we have presented the first full QCD lattice 
calculation involving an arbitrary chemical potential. Furthermore, our results may serve 
as benchmarks for analytic approaches, which can easily tune the quark masses. There 
are many ways forward from here. As a next step, it would be interesting to study the 
cold and dense regime within the current effective theory. Improvements on the present 
state can be made by either carrying the analytic calculations to higher order, including 
additional couplings, or non-perturbatively by inverse Monte Carlo methods along the lines 
of [32, 33]. Finally and most importantly, it is most interesting to explore the possibilities 
of a similar description for the light quark sector, either by extrapolating a higher-order 
hopping expansion or by an alternative formulation within the effective theory context. 
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